stanEx4
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # not see parameters str1..., str2,... using regex as exc='[str1,str2,...]'
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}explanatory variable obs.x have noise
xN(x0.sx),yN(b0+b1*x0,s)
ex8-0.stan
normal regression without explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}ex9.stan
normal regression with explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x10;
}
parameters{
real b0;
real b1;
real<lower=0> s;
real<lower=0> sx;
vector[N] x0;
}
model{
for(i in 1:N){
x[i]~normal(x0[i],sx);
y[i]~normal(b0+b1*x0[i],s);
}
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x0[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] x1;
vector[N1] y1;
for(i in 1:N1){
x1[i]=normal_rng(x10[i],sx);
m1[i]=b0+b1*x10[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)
n1=10
#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -412584
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22)
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22)
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 37 -38.594 0.00664654 0.00175984 1 1 93
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -38.59
## b0 5.57
## b1 2.05
## s 4.18
## m0[1] 7.76
## m0[2] 28.44
## m0[3] 29.48
## m0[4] 11.26
## m0[5] 33.90
## m0[6] 21.13
## m0[7] 20.55
## m0[8] 42.69
## m0[9] 17.27
## m0[10] 30.74
## m0[11] 17.52
## m0[12] 18.34
## m0[13] 31.84
## m0[14] 17.17
## m0[15] 30.77
## m0[16] 19.29
## m0[17] 30.10
## m0[18] 17.60
## m0[19] 12.71
## m0[20] 32.29
## y0[1] 2.17
## y0[2] 27.18
## y0[3] 29.98
## y0[4] 10.35
## y0[5] 31.45
## y0[6] 19.69
## y0[7] 21.51
## y0[8] 41.30
## y0[9] 21.35
## y0[10] 29.66
## y0[11] 12.77
## y0[12] 15.11
## y0[13] 25.14
## y0[14] 17.04
## y0[15] 28.45
## y0[16] 11.80
## y0[17] 34.46
## y0[18] 12.21
## y0[19] 16.23
## y0[20] 29.50
## m1[1] 7.76
## m1[2] 11.65
## m1[3] 15.53
## m1[4] 19.41
## m1[5] 23.29
## m1[6] 27.17
## m1[7] 31.05
## m1[8] 34.93
## m1[9] 38.81
## m1[10] 42.69
## y1[1] 8.05
## y1[2] 11.99
## y1[3] 17.61
## y1[4] 18.23
## y1[5] 23.03
## y1[6] 27.97
## y1[7] 24.03
## y1[8] 37.41
## y1[9] 37.60
## y1[10] 47.32
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -38.80 -38.37 1.49 1.06 -41.80 -37.32 1.01 370 212
## b0 5.47 5.39 2.67 2.42 1.66 9.66 1.02 248 158
## b1 2.06 2.07 0.27 0.24 1.63 2.47 1.02 284 307
## s 4.75 4.60 0.87 0.76 3.61 6.40 1.00 662 726
## m0[1] 7.68 7.62 2.41 2.20 4.23 11.47 1.02 249 161
## m0[2] 28.46 28.48 1.21 1.15 26.50 30.35 1.00 1909 1223
## m0[3] 29.51 29.52 1.29 1.22 27.44 31.50 1.00 1660 1196
## m0[4] 11.19 11.18 2.00 1.86 8.29 14.39 1.02 255 72
## m0[5] 33.95 33.98 1.69 1.58 31.16 36.57 1.00 791 862
## m0[6] 21.12 21.09 1.13 1.08 19.32 22.94 1.00 477 488
## m0[7] 20.53 20.50 1.16 1.09 18.72 22.40 1.01 431 484
## m0[8] 42.79 42.83 2.69 2.43 38.33 46.87 1.01 455 464
## m0[9] 17.24 17.20 1.39 1.27 15.15 19.46 1.02 305 280
## m0[10] 30.78 30.81 1.39 1.30 28.55 32.90 1.00 1322 968
## m0[11] 17.48 17.44 1.37 1.25 15.42 19.66 1.01 310 313
## m0[12] 18.31 18.28 1.31 1.19 16.33 20.39 1.01 331 368
## m0[13] 31.87 31.92 1.48 1.39 29.46 34.19 1.00 1066 893
## m0[14] 17.13 17.10 1.40 1.27 15.04 19.37 1.02 303 290
## m0[15] 30.80 30.84 1.39 1.30 28.57 32.92 1.00 1317 968
## m0[16] 19.26 19.24 1.24 1.14 17.34 21.25 1.01 363 393
## m0[17] 30.13 30.16 1.33 1.25 28.00 32.17 1.00 1497 1165
## m0[18] 17.56 17.51 1.37 1.24 15.51 19.73 1.01 312 326
## m0[19] 12.65 12.63 1.84 1.70 9.98 15.61 1.02 259 87
## m0[20] 32.33 32.37 1.53 1.43 29.82 34.71 1.00 987 866
## y0[1] 7.75 7.76 5.52 5.03 -1.03 16.90 1.00 865 1537
## y0[2] 28.53 28.41 5.10 4.91 20.31 36.89 1.00 2216 2060
## y0[3] 29.55 29.52 4.88 4.66 21.47 37.22 1.00 2120 2024
## y0[4] 11.17 11.19 5.25 5.07 2.69 19.57 1.01 772 1239
## y0[5] 33.91 34.01 5.13 4.71 25.59 42.23 1.00 1794 1878
## y0[6] 21.17 21.07 4.97 4.94 13.15 29.53 1.00 1667 1624
## y0[7] 20.63 20.52 4.84 4.68 12.80 28.78 1.00 1397 1984
## y0[8] 42.94 42.96 5.43 5.26 34.05 51.68 1.00 1203 1372
## y0[9] 17.26 17.20 4.91 4.87 9.43 25.31 1.00 1886 1658
## y0[10] 30.82 30.77 5.14 5.17 22.34 39.01 1.00 1752 1816
## y0[11] 17.46 17.55 5.05 4.92 9.15 25.73 1.00 1586 1591
## y0[12] 18.32 18.35 4.98 4.79 10.12 26.25 1.00 1623 1805
## y0[13] 31.74 31.80 5.09 4.85 23.29 39.61 1.00 1914 1617
## y0[14] 17.04 17.17 5.13 4.97 8.86 25.46 1.00 1595 1728
## y0[15] 30.83 30.97 5.04 4.95 22.58 38.66 1.00 1881 1804
## y0[16] 19.10 19.23 5.00 4.72 10.81 27.10 1.00 1491 2055
## y0[17] 30.19 30.33 5.02 4.90 22.07 38.22 1.00 1983 1783
## y0[18] 17.49 17.50 5.07 4.85 9.23 25.61 1.00 1515 1665
## y0[19] 12.75 12.66 5.30 5.17 4.16 21.31 1.00 953 1475
## y0[20] 32.38 32.36 5.07 4.95 23.94 40.75 1.00 1868 1933
## m1[1] 7.68 7.62 2.41 2.20 4.23 11.47 1.02 249 161
## m1[2] 11.58 11.57 1.96 1.81 8.75 14.71 1.02 256 72
## m1[3] 15.48 15.45 1.55 1.41 13.18 17.96 1.02 279 167
## m1[4] 19.38 19.36 1.23 1.14 17.47 21.36 1.01 368 393
## m1[5] 23.28 23.27 1.07 1.01 21.57 25.00 1.00 808 865
## m1[6] 27.18 27.19 1.14 1.11 25.35 28.98 1.00 1985 1239
## m1[7] 31.08 31.12 1.41 1.32 28.81 33.25 1.00 1241 939
## m1[8] 34.98 35.03 1.79 1.68 32.00 37.76 1.00 707 766
## m1[9] 38.89 38.91 2.23 2.06 35.13 42.31 1.01 532 493
## m1[10] 42.79 42.83 2.69 2.43 38.33 46.87 1.01 455 464
## y1[1] 7.73 7.82 5.38 5.12 -0.99 16.36 1.00 921 937
## y1[2] 11.47 11.39 5.24 4.98 2.87 20.28 1.00 923 1216
## y1[3] 15.24 15.32 5.14 4.92 6.60 23.72 1.00 1058 1339
## y1[4] 19.76 19.85 5.01 4.74 11.46 27.88 1.00 1502 1687
## y1[5] 23.15 23.35 4.86 4.81 15.26 30.80 1.00 1919 1931
## y1[6] 27.24 27.27 4.83 4.53 19.58 35.03 1.00 1811 1842
## y1[7] 31.13 30.98 5.02 4.59 22.65 39.90 1.00 1893 1883
## y1[8] 34.78 34.88 5.17 4.93 26.14 43.30 1.00 1972 1760
## y1[9] 38.88 38.95 5.21 4.94 30.21 47.34 1.01 1129 1359
## y1[10] 42.70 42.74 5.56 5.27 33.46 51.56 1.00 947 1115
#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)
mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 2 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 4 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 2 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 2 finished in 0.4 seconds.
## Chain 1 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 3 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 4 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 1 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
## Chain 4 finished in 0.5 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.6 seconds.seeMCMC(mcmc,'[m,x]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -44.95 -46.14 10.26 8.76 -59.90 -25.84 1.04 118 181
## b0 5.83 5.92 2.58 2.50 1.69 10.02 1.00 749 993
## b1 2.02 2.02 0.27 0.26 1.58 2.46 1.00 675 738
## s 3.58 3.63 1.37 1.38 1.36 5.77 1.03 175 130
## sx 1.46 1.43 0.82 0.89 0.31 2.84 1.05 70 171
## x0[1] 1.02 1.07 1.17 0.92 -1.09 2.91 1.02 1280 803
## x0[2] 12.02 11.83 1.28 1.15 10.23 14.33 1.01 547 922
## x0[3] 13.51 13.17 1.85 1.96 11.20 16.93 1.02 265 639
## x0[4] 3.55 3.42 1.20 1.13 1.78 5.60 1.01 608 2135
## x0[5] 13.38 13.46 1.14 0.97 11.53 15.12 1.01 1097 1684
## x0[6] 6.96 7.11 1.14 1.01 4.96 8.58 1.01 693 1061
## x0[7] 6.30 6.45 1.31 1.26 4.05 8.13 1.02 456 913
## x0[8] 19.05 18.70 1.52 1.22 17.17 22.08 1.01 524 535
## x0[9] 4.92 5.11 1.25 1.10 2.61 6.68 1.01 516 861
## x0[10] 12.10 12.13 1.05 0.86 10.41 13.76 1.01 1459 1482
## x0[11] 5.76 5.81 1.06 0.81 3.96 7.38 1.01 2777 1361
## x0[12] 6.80 6.67 1.12 1.01 5.14 8.77 1.01 851 1908
## x0[13] 13.88 13.62 1.39 1.27 12.08 16.52 1.02 421 870
## x0[14] 5.92 5.88 1.05 0.90 4.25 7.64 1.01 1909 1918
## x0[15] 12.06 12.10 1.09 0.89 10.20 13.90 1.01 2147 1661
## x0[16] 6.42 6.52 1.04 0.84 4.58 8.03 1.01 2140 1140
## x0[17] 10.39 10.60 1.56 1.72 7.70 12.44 1.03 177 1069
## x0[18] 6.80 6.64 1.24 1.23 5.04 8.89 1.01 482 1396
## x0[19] 2.82 3.03 1.25 0.99 0.54 4.57 1.01 736 897
## x0[20] 11.59 11.83 1.47 1.60 9.02 13.54 1.02 268 1117
## m0[1] 8.01 8.03 2.48 2.30 3.89 12.04 1.00 2213 2134
## m0[2] 30.06 30.02 2.46 2.52 26.17 34.13 1.01 613 1630
## m0[3] 33.02 32.74 3.46 4.05 27.88 38.66 1.02 269 1754
## m0[4] 13.08 13.11 2.59 2.64 8.84 17.19 1.01 517 1923
## m0[5] 32.82 32.76 2.40 2.32 29.02 36.89 1.01 1271 2577
## m0[6] 19.95 20.03 2.26 2.22 16.23 23.46 1.01 814 2204
## m0[7] 18.63 18.74 2.57 2.70 14.37 22.62 1.01 474 1729
## m0[8] 44.16 44.35 2.91 2.86 39.12 48.63 1.01 760 2547
## m0[9] 15.85 15.91 2.47 2.52 11.90 19.80 1.01 666 3058
## m0[10] 30.25 30.21 2.15 1.93 26.79 33.94 1.00 1964 2530
## m0[11] 17.53 17.56 2.07 1.84 14.04 20.95 1.00 3322 2484
## m0[12] 19.59 19.58 2.29 2.29 16.00 23.25 1.01 747 2154
## m0[13] 33.78 33.77 2.63 2.82 29.59 38.02 1.01 436 2279
## m0[14] 17.83 17.88 2.15 2.02 14.24 21.31 1.00 1835 2376
## m0[15] 30.18 30.17 2.26 2.11 26.52 33.88 1.00 3026 2206
## m0[16] 18.86 18.87 2.09 1.87 15.37 22.18 1.00 3041 2615
## m0[17] 26.85 27.07 3.22 3.64 21.64 31.71 1.03 165 1437
## m0[18] 19.58 19.50 2.54 2.68 15.54 23.67 1.02 463 2091
## m0[19] 11.63 11.57 2.44 2.40 7.65 15.68 1.00 1064 2328
## m0[20] 29.27 29.37 3.13 3.52 24.23 34.14 1.02 279 1571
## y0[1] 8.06 8.02 4.54 4.00 0.72 15.57 1.00 2867 2468
## y0[2] 30.03 30.43 4.59 4.19 21.95 36.85 1.00 1947 3122
## y0[3] 33.02 33.42 5.19 5.18 24.02 40.59 1.02 444 2752
## y0[4] 13.04 13.38 4.63 4.12 4.90 20.25 1.00 1877 2511
## y0[5] 32.93 32.62 4.44 3.87 25.96 40.65 1.00 2237 3209
## y0[6] 19.96 19.67 4.51 4.06 12.75 27.88 1.00 2369 2697
## y0[7] 18.62 18.14 4.66 4.32 11.78 26.77 1.00 1626 3471
## y0[8] 44.22 44.57 4.79 4.29 36.06 51.48 1.00 2031 2766
## y0[9] 15.87 15.50 4.54 4.10 9.00 23.65 1.00 1788 3538
## y0[10] 30.25 30.24 4.39 3.90 22.99 37.41 1.00 3540 3092
## y0[11] 17.57 17.55 4.34 3.86 10.31 24.83 1.00 3680 3348
## y0[12] 19.55 19.80 4.49 3.95 11.80 26.47 1.00 2457 3358
## y0[13] 33.90 34.31 4.54 4.23 25.78 40.50 1.00 1507 2918
## y0[14] 17.85 18.02 4.37 3.83 10.47 24.71 1.00 3754 3301
## y0[15] 30.13 29.98 4.34 3.83 23.12 37.43 1.00 4201 2938
## y0[16] 18.75 18.65 4.36 3.85 11.64 25.97 1.00 4073 3483
## y0[17] 26.75 26.22 4.95 4.83 19.62 35.48 1.01 692 1964
## y0[18] 19.58 20.03 4.61 4.37 11.63 26.47 1.00 1569 2643
## y0[19] 11.61 11.30 4.54 4.14 4.55 19.26 1.00 2653 3083
## y0[20] 29.23 28.78 4.96 4.86 22.14 38.09 1.01 665 2501
## m1[1] 8.00 8.05 2.33 2.25 4.30 11.77 1.00 771 1011
## m1[2] 11.82 11.84 1.90 1.81 8.77 14.93 1.00 838 1231
## m1[3] 15.64 15.67 1.53 1.45 13.17 18.15 1.00 994 1394
## m1[4] 19.47 19.46 1.25 1.17 17.44 21.49 1.00 1354 1500
## m1[5] 23.29 23.30 1.13 1.08 21.44 25.14 1.00 1797 1822
## m1[6] 27.11 27.12 1.24 1.19 25.08 29.15 1.00 1432 2294
## m1[7] 30.94 30.96 1.51 1.44 28.44 33.41 1.00 1060 1840
## m1[8] 34.76 34.76 1.88 1.81 31.66 37.84 1.00 891 1468
## m1[9] 38.58 38.56 2.31 2.23 34.80 42.34 1.00 806 1186
## m1[10] 42.41 42.39 2.76 2.68 37.88 46.87 1.00 763 1018
## x1[1] 1.07 1.09 1.67 1.18 -1.66 3.73 1.02 3843 2173
## x1[2] 2.97 2.96 1.69 1.14 0.20 5.81 1.02 3750 1971
## x1[3] 4.87 4.85 1.68 1.15 2.17 7.64 1.01 3855 2675
## x1[4] 6.73 6.73 1.70 1.18 4.06 9.41 1.02 3709 2192
## x1[5] 8.63 8.62 1.65 1.14 5.90 11.28 1.02 4148 2087
## x1[6] 10.53 10.53 1.66 1.14 7.77 13.28 1.02 4013 2495
## x1[7] 12.46 12.45 1.68 1.18 9.80 15.16 1.02 4024 1813
## x1[8] 14.35 14.34 1.66 1.13 11.59 16.95 1.02 4028 2139
## x1[9] 16.23 16.20 1.68 1.18 13.57 19.08 1.02 3892 2389
## x1[10] 18.11 18.12 1.67 1.14 15.46 20.83 1.02 3962 2205
## y1[1] 7.95 7.80 4.59 4.13 0.62 15.38 1.00 1619 2319
## y1[2] 11.82 11.78 4.34 3.88 4.60 18.92 1.00 2270 1950
## y1[3] 15.68 15.71 4.16 3.60 9.03 22.43 1.00 3144 3702
## y1[4] 19.45 19.45 3.88 3.40 12.96 25.86 1.00 3779 3648
## y1[5] 23.20 23.28 3.99 3.43 16.55 29.73 1.00 3594 2668
## y1[6] 26.99 27.03 4.11 3.63 20.31 33.64 1.00 3395 3368
## y1[7] 30.93 30.91 4.10 3.61 24.11 37.57 1.00 3166 2843
## y1[8] 34.76 34.80 4.23 3.78 27.89 41.60 1.00 2548 2759
## y1[9] 38.75 38.71 4.49 4.10 31.52 45.92 1.00 1714 2471
## y1[10] 42.44 42.44 4.60 4.26 34.90 50.07 1.00 1454 2259
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
outlier
ex10.stan
objective variable have outlier, y~cauchy(b0+b1*x,s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~cauchy(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=cauchy_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=cauchy_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -8140.49
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 20 -32.945 0.000939809 2.02029e-05 1 1 39
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -32.94
## b0 5.88
## b1 1.97
## s 3.15
## m0[1] 11.78
## m0[2] 10.27
## m0[3] 20.62
## m0[4] 21.91
## m0[5] 13.05
## m0[6] 6.70
## m0[7] 10.94
## m0[8] 18.87
## m0[9] 12.98
## m0[10] 20.66
## m0[11] 22.87
## m0[12] 6.77
## m0[13] 13.26
## m0[14] 16.58
## m0[15] 23.05
## m0[16] 13.73
## m0[17] 12.55
## m0[18] 6.07
## m0[19] 17.99
## m0[20] 11.58
## y0[1] 9.84
## y0[2] 9.29
## y0[3] 18.63
## y0[4] 18.17
## y0[5] 13.89
## y0[6] 8.13
## y0[7] 16.62
## y0[8] 14.74
## y0[9] 11.95
## y0[10] 23.93
## y0[11] 26.61
## y0[12] 4.81
## y0[13] 10.59
## y0[14] 16.98
## y0[15] 32.01
## y0[16] 12.84
## y0[17] 11.85
## y0[18] 7.09
## y0[19] 8.88
## y0[20] 10.60
## m1[1] 6.07
## m1[2] 7.95
## m1[3] 9.84
## m1[4] 11.73
## m1[5] 13.61
## m1[6] 15.50
## m1[7] 17.39
## m1[8] 19.27
## m1[9] 21.16
## m1[10] 23.05
## y1[1] 12.71
## y1[2] 5.00
## y1[3] 6.60
## y1[4] 15.00
## y1[5] 13.39
## y1[6] 15.26
## y1[7] 17.53
## y1[8] 20.62
## y1[9] 24.12
## y1[10] 25.35
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -33.47 -33.16 1.34 1.16 -35.96 -31.96 1.01 529 714
## b0 5.81 5.80 1.67 1.57 2.95 8.57 1.01 299 283
## b1 1.98 1.99 0.31 0.30 1.46 2.51 1.01 393 604
## s 3.57 3.49 0.65 0.60 2.73 4.72 1.00 1470 1506
## m0[1] 11.75 11.72 0.98 0.96 10.17 13.40 1.01 387 477
## m0[2] 10.23 10.22 1.12 1.07 8.43 12.12 1.01 334 322
## m0[3] 20.65 20.66 1.24 1.19 18.61 22.68 1.00 1549 1383
## m0[4] 21.95 21.96 1.40 1.34 19.67 24.25 1.00 1225 1308
## m0[5] 13.03 13.01 0.89 0.87 11.58 14.52 1.00 484 604
## m0[6] 6.63 6.63 1.56 1.48 4.01 9.24 1.01 300 276
## m0[7] 10.90 10.88 1.06 1.01 9.19 12.69 1.01 351 360
## m0[8] 18.89 18.88 1.05 1.00 17.12 20.58 1.00 1905 1448
## m0[9] 12.95 12.93 0.90 0.88 11.50 14.45 1.00 475 600
## m0[10] 20.69 20.70 1.24 1.19 18.64 22.72 1.00 1541 1383
## m0[11] 22.91 22.94 1.52 1.45 20.43 25.40 1.00 1027 1290
## m0[12] 6.70 6.71 1.55 1.47 4.10 9.29 1.01 300 282
## m0[13] 13.24 13.22 0.88 0.86 11.82 14.71 1.00 511 635
## m0[14] 16.59 16.59 0.88 0.88 15.10 18.02 1.00 1760 1513
## m0[15] 23.09 23.12 1.55 1.47 20.57 25.61 1.00 995 1290
## m0[16] 13.72 13.70 0.86 0.86 12.34 15.14 1.00 580 649
## m0[17] 12.53 12.50 0.92 0.90 11.02 14.08 1.00 436 534
## m0[18] 5.99 5.99 1.65 1.55 3.19 8.73 1.01 299 274
## m0[19] 18.01 17.99 0.97 0.95 16.38 19.58 1.00 1957 1542
## m0[20] 11.55 11.52 1.00 0.97 9.94 13.23 1.01 377 461
## y0[1] 11.77 11.77 3.80 3.62 5.46 18.03 1.00 1489 1826
## y0[2] 10.20 10.19 3.81 3.62 4.07 16.37 1.00 1256 1692
## y0[3] 20.53 20.60 3.83 3.77 14.22 26.55 1.00 1954 2012
## y0[4] 21.98 22.05 3.91 3.85 15.63 28.13 1.00 1746 1972
## y0[5] 13.01 13.03 3.72 3.63 7.06 19.00 1.00 1752 1794
## y0[6] 6.59 6.62 4.11 3.96 -0.48 13.03 1.00 947 1720
## y0[7] 10.99 11.01 3.75 3.62 4.84 17.18 1.00 1934 1382
## y0[8] 18.87 18.96 3.76 3.70 12.59 25.24 1.00 2068 1905
## y0[9] 12.94 12.99 3.73 3.58 6.61 18.96 1.00 1930 1997
## y0[10] 20.71 20.69 3.94 3.66 14.12 27.19 1.00 1845 1675
## y0[11] 22.85 22.83 4.08 3.75 16.36 29.65 1.00 1691 1483
## y0[12] 6.67 6.82 3.94 3.81 0.18 13.14 1.00 1201 1475
## y0[13] 13.24 13.29 3.67 3.56 7.23 19.13 1.00 1914 1968
## y0[14] 16.62 16.57 3.72 3.67 10.53 22.54 1.00 1799 1893
## y0[15] 23.04 22.99 4.01 3.92 16.53 29.69 1.00 1943 1893
## y0[16] 13.54 13.54 3.71 3.53 7.58 19.58 1.00 1914 1732
## y0[17] 12.49 12.50 3.75 3.68 6.37 18.78 1.00 1713 1587
## y0[18] 5.89 6.02 4.08 3.76 -1.03 12.10 1.00 1038 1407
## y0[19] 18.05 18.03 3.77 3.58 11.97 24.11 1.00 2095 1972
## y0[20] 11.47 11.54 3.81 3.66 5.23 17.59 1.00 1556 1650
## m1[1] 5.99 5.99 1.65 1.55 3.19 8.73 1.01 299 274
## m1[2] 7.89 7.91 1.40 1.33 5.60 10.25 1.01 306 275
## m1[3] 9.79 9.79 1.17 1.12 7.94 11.74 1.01 326 314
## m1[4] 11.69 11.66 0.98 0.97 10.10 13.36 1.01 384 477
## m1[5] 13.59 13.58 0.87 0.86 12.20 15.03 1.00 561 634
## m1[6] 15.49 15.50 0.85 0.83 14.14 16.87 1.00 1091 1140
## m1[7] 17.39 17.39 0.93 0.92 15.81 18.90 1.00 1916 1554
## m1[8] 19.29 19.28 1.09 1.03 17.45 21.06 1.00 1869 1518
## m1[9] 21.19 21.20 1.30 1.25 19.06 23.34 1.00 1440 1364
## m1[10] 23.09 23.12 1.55 1.47 20.57 25.61 1.00 995 1290
## y1[1] 5.87 5.81 4.00 3.88 -0.54 12.36 1.00 1204 1512
## y1[2] 7.95 7.94 3.86 3.85 1.88 14.40 1.00 970 1846
## y1[3] 9.78 9.82 3.76 3.63 3.48 15.98 1.00 1561 1755
## y1[4] 11.63 11.70 3.83 3.80 5.34 17.81 1.00 1763 1658
## y1[5] 13.70 13.67 3.78 3.59 7.69 19.94 1.00 1812 1892
## y1[6] 15.39 15.35 3.67 3.48 9.52 21.54 1.00 1891 1931
## y1[7] 17.27 17.30 3.64 3.52 11.35 23.24 1.00 1797 1880
## y1[8] 19.30 19.20 3.65 3.49 13.46 25.33 1.00 2026 2013
## y1[9] 21.25 21.32 3.92 3.81 14.84 27.58 1.00 1943 1822
## y1[10] 23.20 23.23 3.98 3.70 16.82 29.83 1.00 2004 1839
mdl=cmdstan_model('./ex10.stan')
fn(mdl,data)## Initial log joint probability = -124.542
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 24 -9.48434 0.000122605 0.00219896 0.9709 0.9709 34
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -9.48
## b0 4.45
## b1 2.07
## s 0.47
## m0[1] 10.66
## m0[2] 9.07
## m0[3] 19.95
## m0[4] 21.30
## m0[5] 11.99
## m0[6] 5.31
## m0[7] 9.77
## m0[8] 18.11
## m0[9] 11.92
## m0[10] 19.99
## m0[11] 22.31
## m0[12] 5.39
## m0[13] 12.22
## m0[14] 15.71
## m0[15] 22.50
## m0[16] 12.71
## m0[17] 11.47
## m0[18] 4.65
## m0[19] 17.19
## m0[20] 10.45
## y0[1] 11.60
## y0[2] 8.74
## y0[3] 21.60
## y0[4] 21.65
## y0[5] 11.81
## y0[6] 8.40
## y0[7] 9.61
## y0[8] 18.20
## y0[9] 12.06
## y0[10] 19.18
## y0[11] 19.53
## y0[12] 4.47
## y0[13] 12.48
## y0[14] 15.24
## y0[15] 21.71
## y0[16] 12.76
## y0[17] 11.65
## y0[18] 4.39
## y0[19] 19.67
## y0[20] 9.57
## m1[1] 4.65
## m1[2] 6.63
## m1[3] 8.62
## m1[4] 10.60
## m1[5] 12.58
## m1[6] 14.57
## m1[7] 16.55
## m1[8] 18.54
## m1[9] 20.52
## m1[10] 22.50
## y1[1] 4.61
## y1[2] 7.05
## y1[3] 8.36
## y1[4] 12.84
## y1[5] 12.42
## y1[6] 15.57
## y1[7] 17.02
## y1[8] 25.72
## y1[9] 19.62
## y1[10] 20.69
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -11.99 -11.61 1.44 1.16 -15.02 -10.42 1.01 531 594
## b0 4.45 4.45 0.37 0.33 3.84 5.07 1.00 669 626
## b1 2.08 2.08 0.07 0.06 1.98 2.20 1.00 777 1001
## s 0.60 0.58 0.20 0.19 0.34 0.97 1.01 525 286
## m0[1] 10.70 10.69 0.24 0.23 10.33 11.11 1.01 853 766
## m0[2] 9.11 9.10 0.26 0.25 8.69 9.55 1.00 759 685
## m0[3] 20.07 20.03 0.30 0.27 19.64 20.63 1.00 1702 1380
## m0[4] 21.43 21.39 0.33 0.29 20.96 22.05 1.00 1540 1299
## m0[5] 12.05 12.03 0.22 0.22 11.71 12.45 1.01 995 1003
## m0[6] 5.32 5.31 0.35 0.32 4.74 5.90 1.00 675 621
## m0[7] 9.81 9.80 0.25 0.24 9.42 10.24 1.00 794 699
## m0[8] 18.21 18.19 0.26 0.23 17.84 18.70 1.00 1764 1378
## m0[9] 11.97 11.96 0.22 0.22 11.63 12.37 1.01 984 1003
## m0[10] 20.11 20.07 0.30 0.27 19.68 20.67 1.00 1696 1380
## m0[11] 22.45 22.41 0.36 0.32 21.94 23.11 1.00 1440 1253
## m0[12] 5.40 5.39 0.34 0.32 4.82 5.97 1.00 675 621
## m0[13] 12.28 12.26 0.22 0.21 11.94 12.67 1.01 1026 1047
## m0[14] 15.79 15.77 0.23 0.21 15.47 16.20 1.00 1656 1287
## m0[15] 22.64 22.60 0.36 0.32 22.12 23.31 1.00 1422 1258
## m0[16] 12.77 12.75 0.22 0.21 12.44 13.16 1.01 1103 1207
## m0[17] 11.52 11.51 0.23 0.22 11.17 11.93 1.01 932 896
## m0[18] 4.65 4.65 0.36 0.33 4.05 5.26 1.00 670 626
## m0[19] 17.29 17.26 0.25 0.22 16.94 17.74 1.00 1755 1436
## m0[20] 10.49 10.48 0.24 0.23 10.11 10.91 1.01 837 709
## y0[1] 10.60 10.68 20.96 0.90 6.13 14.42 1.00 1809 2015
## y0[2] 9.37 9.10 9.44 0.98 5.27 13.39 1.00 1915 1838
## y0[3] 19.68 20.05 31.14 1.01 15.90 24.55 1.00 1971 1892
## y0[4] 26.33 21.44 193.08 0.93 18.01 25.25 1.00 1930 1793
## y0[5] 12.50 12.03 19.59 0.91 8.72 16.22 1.00 2048 1931
## y0[6] 5.46 5.31 10.35 1.03 1.64 9.50 1.00 1704 1787
## y0[7] 9.28 9.82 31.06 0.95 5.65 13.58 1.00 1936 1843
## y0[8] 18.69 18.16 16.42 0.90 14.64 22.03 1.00 2067 1877
## y0[9] 12.28 11.94 20.98 0.94 8.13 15.61 1.00 1795 1973
## y0[10] 19.98 20.07 17.56 0.97 15.56 23.94 1.00 1901 1974
## y0[11] 22.25 22.40 8.54 0.98 18.75 26.16 1.00 1931 1857
## y0[12] 2.80 5.42 79.01 0.98 1.74 9.30 1.00 1868 2052
## y0[13] 12.18 12.27 9.12 0.96 8.38 15.68 1.00 1980 1980
## y0[14] 15.82 15.77 10.70 0.85 11.99 19.71 1.00 2129 2036
## y0[15] 23.04 22.60 18.71 0.95 18.65 26.72 1.00 1975 1822
## y0[16] 12.75 12.79 45.54 0.87 8.75 16.04 1.00 1995 1694
## y0[17] 11.91 11.54 23.89 0.92 7.96 15.45 1.00 2089 1835
## y0[18] 4.36 4.65 23.74 1.00 0.73 8.31 1.00 1622 1755
## y0[19] 16.60 17.25 36.10 0.85 13.46 21.05 1.00 2140 1972
## y0[20] 10.63 10.49 7.08 0.90 7.49 14.11 1.00 1927 1787
## m1[1] 4.65 4.65 0.36 0.33 4.05 5.26 1.00 670 626
## m1[2] 6.65 6.64 0.31 0.29 6.13 7.17 1.00 691 612
## m1[3] 8.65 8.64 0.27 0.26 8.22 9.10 1.00 740 685
## m1[4] 10.65 10.63 0.24 0.23 10.27 11.06 1.01 849 766
## m1[5] 12.65 12.63 0.22 0.21 12.31 13.03 1.01 1081 1137
## m1[6] 14.64 14.62 0.22 0.21 14.32 15.03 1.00 1458 1342
## m1[7] 16.64 16.61 0.24 0.22 16.31 17.08 1.00 1721 1323
## m1[8] 18.64 18.61 0.27 0.24 18.26 19.14 1.00 1762 1319
## m1[9] 20.64 20.60 0.31 0.28 20.20 21.22 1.00 1630 1322
## m1[10] 22.64 22.60 0.36 0.32 22.12 23.31 1.00 1422 1258
## y1[1] 5.24 4.61 28.58 0.94 0.41 7.86 1.00 1578 1631
## y1[2] 6.35 6.64 19.10 0.91 2.88 10.33 1.00 1602 1812
## y1[3] 9.78 8.64 24.87 0.97 5.15 12.72 1.00 1801 1987
## y1[4] 10.90 10.66 34.15 0.94 6.65 14.53 1.00 1942 1933
## y1[5] 12.92 12.63 19.90 0.90 8.44 16.48 1.00 1945 1894
## y1[6] 15.50 14.65 49.43 0.93 10.52 18.58 1.00 1947 1896
## y1[7] 16.40 16.64 16.06 0.88 12.98 20.80 1.00 2082 1707
## y1[8] 17.40 18.62 29.80 0.91 14.62 22.55 1.00 1867 1784
## y1[9] 20.66 20.62 12.96 0.95 16.71 24.68 1.00 1924 1969
## y1[10] 17.45 22.62 193.75 1.00 19.09 26.13 1.00 2097 1755
censored data
objective variable has NA
(x,y) i=1-N
(x0,y0) i=1-N0
x1 i=1-N1, y1=NA
(x,y)~N((mx,my),(sx2,sy2,sxy))
(x0,y0)~N((mx,my),(sx2,sy2,sxy))
x1~N(mx,sx2)ex11-0.stan
data{
int N0;
array[N0] vector[2] xy;
int N1;
vector[N1] x1;
}
parameters{
vector[2] m;
cov_matrix[2] s;
}
model{
target+=multi_normal_lpdf(xy | m, s);
x1~normal(m[1],s[1,1]^.5);
}
generated quantities{
vector[2] xy1;
xy1=multi_normal_rng(m,s);
real cr;
cr=s[1,2]/(s[1,1]*s[2,2])^.5;
}n=30
x=runif(n,0,9)
y=rnorm(n,10+3*x,4)
cor(x,y)## [1] 0.905qplot(x,y)
L=4
n0=sum(x>L)
x0=x[x>L]
y0=y[x>L]
x1=x[x<=L]
n1=sum(x<=L)
cor(x0,y0)## [1] 0.716qplot(x0,y0)
mdl=cmdstan_model('./ex11-0.stan')
data=list(N0=n0,xy=matrix(c(x0,y0),ncol=2),N1=n1,x1=x1)
mle=mdl$optimize(data=data)## Initial log joint probability = -107248
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 59 -106.626 0.000515973 0.0018835 1 1 79
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -106.63
## m[1] 4.61
## m[2] 25.09
## s[1,1] 7.31
## s[2,1] 20.79
## s[1,2] 20.79
## s[2,2] 72.51
## xy1[1] 5.82
## xy1[2] 25.84
## cr 0.90mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 finished in 0.4 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.5 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -102.16 -101.85 1.62 1.46 -105.22 -100.13 1.00 664 1041
## m[1] 4.61 4.62 0.54 0.50 3.68 5.49 1.00 726 708
## m[2] 25.16 25.26 2.56 2.35 20.90 29.31 1.01 400 374
## s[1,1] 9.43 8.92 2.83 2.50 5.80 14.79 1.01 1034 1347
## s[2,1] 26.35 24.53 10.92 9.61 11.86 46.80 1.01 437 608
## s[1,2] 26.35 24.53 10.92 9.61 11.86 46.80 1.01 437 608
## s[2,2] 100.27 88.62 49.64 40.97 41.46 198.19 1.01 425 517
## xy1[1] 4.60 4.61 3.14 2.99 -0.61 9.70 1.00 1631 1673
## xy1[2] 25.27 25.60 10.35 9.43 8.40 41.68 1.00 1468 1346
## cr 0.86 0.89 0.12 0.07 0.66 0.96 1.01 529 630
xy=mcmc$draws('xy1')
cor(xy[,,1],xy[,,2])## [1] 0.835qplot(xy[,,1],xy[,,2])
objective variable is censored
y i=1-N, y~N(m,s)
actual ya i=1-Na
lower censored yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
upper censored yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)
cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))
P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)ex11-1.stan
data{
int N;
vector[N] x;
vector[N] y;
real L;
int Nl;
vector[Nl] xl;
real U;
int Nu;
vector[Nu] xu;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
for(i in 1:Nl)
target+=normal_lcdf(L | b0+b1*xl[i],s);
for(i in 1:Nu)
target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)
xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitionmdl=cmdstan_model('./ex8-0.stan')
data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mle=mdl$optimize(data=data)## Initial log joint probability = -19025.6
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## 24 -7.43161 0.000242549 0.00045938 0.9732 0.9732 53
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -7.43
## b0 15.10
## b1 1.94
## s 1.39
## m0[1] 20.06
## m0[2] 26.77
## m0[3] 28.06
## m0[4] 24.88
## m0[5] 18.24
## m0[6] 19.16
## m0[7] 28.15
## m0[8] 20.52
## m0[9] 19.25
## y0[1] 20.67
## y0[2] 26.52
## y0[3] 30.90
## y0[4] 26.40
## y0[5] 19.86
## y0[6] 20.52
## y0[7] 27.71
## y0[8] 20.59
## y0[9] 18.69
## m1[1] 15.79
## m1[2] 17.61
## m1[3] 19.42
## m1[4] 21.24
## m1[5] 23.05
## m1[6] 24.87
## m1[7] 26.69
## m1[8] 28.50
## m1[9] 30.32
## m1[10] 32.13
## y1[1] 15.53
## y1[2] 17.54
## y1[3] 20.63
## y1[4] 21.65
## y1[5] 24.75
## y1[6] 26.84
## y1[7] 26.62
## y1[8] 29.61
## y1[9] 29.45
## y1[10] 32.79mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -8.84 -8.48 1.48 1.19 -11.78 -7.27 1.00 446 310
## b0 15.17 15.16 1.48 1.32 12.89 17.49 1.02 270 315
## b1 1.93 1.93 0.33 0.29 1.40 2.46 1.02 342 463
## s 1.93 1.79 0.68 0.53 1.19 3.14 1.01 605 526
## m0[1] 20.11 20.09 0.83 0.70 18.86 21.47 1.02 335 311
## m0[2] 26.78 26.77 0.98 0.87 25.30 28.35 1.00 1491 789
## m0[3] 28.06 28.05 1.14 1.01 26.30 29.89 1.00 1126 684
## m0[4] 24.90 24.90 0.79 0.67 23.72 26.16 1.00 1859 759
## m0[5] 18.29 18.27 1.04 0.90 16.70 19.91 1.02 279 327
## m0[6] 19.21 19.19 0.93 0.80 17.81 20.65 1.02 292 333
## m0[7] 28.16 28.14 1.16 1.02 26.37 30.00 1.00 1106 685
## m0[8] 20.57 20.54 0.79 0.66 19.37 21.86 1.01 382 388
## m0[9] 19.30 19.27 0.92 0.79 17.91 20.72 1.02 293 335
## y0[1] 20.19 20.15 2.19 1.83 16.73 23.70 1.00 1677 1222
## y0[2] 26.79 26.76 2.31 1.96 23.11 30.36 1.00 1562 1234
## y0[3] 28.04 28.06 2.28 2.05 24.30 31.68 1.00 1718 1417
## y0[4] 24.87 24.81 2.21 1.90 21.45 28.32 1.00 1955 1226
## y0[5] 18.21 18.24 2.29 1.99 14.61 21.75 1.00 1145 1296
## y0[6] 19.30 19.26 2.32 1.91 15.71 22.75 1.00 960 1186
## y0[7] 28.11 28.15 2.36 2.05 24.08 31.79 1.00 1761 1595
## y0[8] 20.59 20.54 2.26 1.79 17.26 23.97 1.00 1279 1188
## y0[9] 19.20 19.20 2.19 1.94 15.63 22.61 1.00 1122 1257
## m1[1] 15.86 15.84 1.38 1.21 13.76 18.02 1.02 270 316
## m1[2] 17.67 17.64 1.12 0.97 15.93 19.40 1.02 274 321
## m1[3] 19.47 19.45 0.90 0.78 18.11 20.88 1.02 297 335
## m1[4] 21.28 21.22 0.74 0.60 20.18 22.51 1.01 461 382
## m1[5] 23.08 23.06 0.70 0.58 22.06 24.20 1.00 1072 552
## m1[6] 24.89 24.88 0.78 0.67 23.71 26.15 1.00 1860 759
## m1[7] 26.70 26.69 0.97 0.86 25.23 28.24 1.01 1510 795
## m1[8] 28.50 28.49 1.20 1.06 26.63 30.42 1.01 1028 703
## m1[9] 30.31 30.29 1.47 1.26 27.99 32.65 1.01 788 580
## m1[10] 32.11 32.09 1.75 1.51 29.39 34.90 1.01 668 636
## y1[1] 15.77 15.71 2.48 2.14 11.80 19.51 1.01 718 934
## y1[2] 17.63 17.62 2.34 2.01 13.98 21.31 1.01 990 1194
## y1[3] 19.54 19.46 2.28 1.97 15.97 23.22 1.00 1292 1602
## y1[4] 21.33 21.32 2.19 1.92 17.98 24.74 1.00 1697 1700
## y1[5] 23.10 23.10 2.14 1.89 19.73 26.63 1.00 1839 1656
## y1[6] 24.83 24.82 2.27 1.89 21.45 28.08 1.00 1758 1817
## y1[7] 26.61 26.58 2.32 1.96 22.98 30.21 1.00 1839 1662
## y1[8] 28.52 28.49 2.46 2.06 24.69 32.42 1.00 1790 1277
## y1[9] 30.30 30.35 2.42 2.09 26.47 34.01 1.00 1889 1485
## y1[10] 32.17 32.20 2.84 2.40 27.89 36.41 1.00 1161 946
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
data=list(N=n,x=xya$x,y=xya$y,
L=L,Nl=nl,xl=xyl$x,
U=U,Nu=nu,xu=xyu$x,
N1=n1,x1=x1)
mdl=cmdstan_model('./ex11-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -205.788
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 31 -16.1894 8.93464e-05 0.0002329 0.88 0.88 42
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -16.19
## b0 11.35
## b1 2.74
## s 2.27
## m0[1] 18.37
## m0[2] 27.85
## m0[3] 29.68
## m0[4] 25.19
## m0[5] 15.79
## m0[6] 17.09
## m0[7] 29.81
## m0[8] 19.02
## m0[9] 17.22
## y0[1] 17.99
## y0[2] 27.24
## y0[3] 28.81
## y0[4] 26.97
## y0[5] 14.17
## y0[6] 19.78
## y0[7] 29.48
## y0[8] 20.25
## y0[9] 19.30
## m1[1] 12.33
## m1[2] 14.89
## m1[3] 17.46
## m1[4] 20.03
## m1[5] 22.60
## m1[6] 25.17
## m1[7] 27.73
## m1[8] 30.30
## m1[9] 32.87
## m1[10] 35.44
## y1[1] 9.00
## y1[2] 10.57
## y1[3] 21.79
## y1[4] 21.48
## y1[5] 22.23
## y1[6] 24.70
## y1[7] 28.68
## y1[8] 31.69
## y1[9] 32.46
## y1[10] 37.90mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m',ch=T)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -17.19 -16.74 1.71 1.15 -20.16 -15.53 1.02 218 156
## b0 9.94 10.43 2.89 1.70 5.98 12.92 1.03 132 113
## b1 3.10 2.98 0.71 0.40 2.42 4.05 1.02 159 163
## s 3.31 2.97 1.66 0.85 1.90 5.60 1.02 229 208
## m0[1] 17.87 18.08 1.37 0.97 15.69 19.52 1.03 173 151
## m0[2] 28.59 28.33 1.92 1.21 26.60 31.20 1.01 810 233
## m0[3] 30.65 30.30 2.33 1.43 28.38 33.69 1.01 482 216
## m0[4] 25.58 25.45 1.40 0.93 23.93 27.45 1.00 1409 489
## m0[5] 14.96 15.27 1.87 1.21 12.11 16.97 1.03 140 138
## m0[6] 16.43 16.70 1.60 1.09 13.92 18.23 1.03 145 143
## m0[7] 30.80 30.44 2.36 1.45 28.51 33.89 1.01 469 216
## m0[8] 18.61 18.77 1.27 0.93 16.65 20.18 1.03 207 201
## m0[9] 16.57 16.84 1.58 1.08 14.10 18.34 1.03 145 144
## y0[1] 17.84 18.10 3.96 3.11 11.90 23.21 1.00 1152 595
## y0[2] 28.63 28.38 4.13 3.29 23.00 35.11 1.01 1714 814
## y0[3] 30.59 30.20 4.50 3.43 24.60 37.29 1.01 1581 566
## y0[4] 25.56 25.41 3.78 3.25 19.67 31.65 1.00 1830 852
## y0[5] 14.93 15.24 4.33 3.21 8.51 20.65 1.01 754 278
## y0[6] 16.45 16.63 3.92 3.14 10.30 21.76 1.01 1117 279
## y0[7] 30.77 30.51 4.24 3.33 24.73 37.45 1.00 1620 538
## y0[8] 18.69 18.70 4.03 3.07 12.60 24.43 1.00 1787 780
## y0[9] 16.65 16.87 4.19 3.08 10.45 22.55 1.01 905 377
## m1[1] 11.04 11.50 2.66 1.60 7.39 13.81 1.03 132 122
## m1[2] 13.95 14.29 2.06 1.32 10.90 16.14 1.03 137 140
## m1[3] 16.85 17.10 1.53 1.05 14.46 18.59 1.03 146 146
## m1[4] 19.75 19.87 1.15 0.86 17.94 21.26 1.02 363 210
## m1[5] 22.65 22.65 1.09 0.82 21.10 24.17 1.01 1103 712
## m1[6] 25.55 25.43 1.39 0.93 23.91 27.43 1.00 1415 500
## m1[7] 28.45 28.20 1.90 1.20 26.48 31.03 1.01 837 233
## m1[8] 31.36 30.97 2.48 1.51 28.97 34.55 1.01 427 217
## m1[9] 34.26 33.75 3.10 1.85 31.32 38.19 1.01 319 205
## m1[10] 37.16 36.53 3.73 2.18 33.67 42.04 1.01 271 182
## y1[1] 11.04 11.46 4.46 3.52 4.29 16.96 1.01 629 371
## y1[2] 13.86 14.13 4.43 3.39 6.95 20.00 1.01 762 282
## y1[3] 16.86 17.10 3.81 3.09 10.82 22.36 1.01 1111 708
## y1[4] 19.85 19.94 3.84 3.08 14.00 25.46 1.00 1259 856
## y1[5] 22.64 22.60 3.76 3.08 17.20 28.63 1.00 1718 1324
## y1[6] 25.59 25.44 4.12 3.10 19.93 31.55 1.00 2058 1242
## y1[7] 28.52 28.26 4.22 3.27 22.76 34.84 1.01 1542 738
## y1[8] 31.36 31.03 4.39 3.31 25.85 38.01 1.00 1282 398
## y1[9] 34.18 33.69 5.18 3.49 28.10 40.97 1.00 828 344
## y1[10] 37.01 36.35 5.31 3.64 30.96 44.13 1.01 597 390
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
sensitivity/specificity analysis
ex14.stan
estimating sens and spec
data {
int N;
array[N] int x;
array[N] int y;
}
parameters {
real<lower=0,upper=1> p;
real<lower=0,upper=1> se;
real<lower=0,upper=1> sp;
}
model {
p~uniform(0,1);
se~uniform(0,1);
sp~uniform(0,1);
for (i in 1:N) {
y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
}
}
generated quantities {
real ppv;
real npv;
ppv=se*p/((se*p)+((1-p)*(1-sp)));
npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}n=20
data=list(N=n,
x=sample(0:1,n,replace=T),
y=sample(0:1,n,replace=T))
mdl=cmdstan_model('./ex14.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -17.2496
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -13.393 0.000174194 1.81971e-05 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.39
## p 0.72
## se 0.56
## sp 0.64
## ppv 0.79
## npv 0.36mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -19.36 -19.00 1.39 1.17 -22.01 -17.83 1.00 776 876
## p 0.49 0.48 0.29 0.37 0.05 0.95 1.01 693 566
## se 0.54 0.54 0.14 0.15 0.30 0.78 1.00 2394 1653
## sp 0.61 0.62 0.14 0.15 0.38 0.82 1.00 2947 1419
## ppv 0.55 0.57 0.29 0.38 0.06 0.97 1.01 724 627
## npv 0.56 0.59 0.29 0.37 0.06 0.97 1.01 726 548
2x2 cross table
Effect occur y=1 with probabilty p01, p11 from each cause x{0,1}
event frequncy nxy of effect y{0,1} by cause x{0,1}
n01~B(n0.,p0)
n11~B(n1.,p1)
n01=n0p0, n00=n0(1-p0)
n11=n1p1, n10=n1(1-p1)
p00=n00/n=n0(1-p0)/n, p01=n01/n=n0p0/n
p10=n10/n=n1(1-p1)/n, p11=n11/n=n1p1/n
Cramer'V (chi2/n/(min(row,column)-1))^.5
in 2x2
crv =(n11*n00-n10*n01)/(n0.*n1.*n.0*n.1)^.5
=(n0(1-p0)n1p1-n0p0n1(1-p1))/(n0n1(n0(1-p0)+n1(1-p1))(n0p0n1p1))^.5
kappa coefficient k=(po-pe)/(1-pe)
po: Observed agreement (proportion of times both raters agreed)
pe: Expected agreement under independence
po=p00+p11
=(n0(1-p0)+n1p1)/n
pe=(p00+p01)(p00+p10)(p10+p11)(p01+p11)
=n0/n*(n0(1-p0)+n1(1-p1))/n*(n0p0+n1p1)/n*n1/nex16-1.stan
data {
int n;
int n0;
int n01;
int n1;
int n11;
}
parameters {
real<lower=0,upper=1> p0;
real<lower=0,upper=1> p1;
}
model {
n01~binomial(n0,p0);
n11~binomial(n1,p1);
}
generated quantities {
real RR;
RR=p1/p0;
real OR;
OR=(p1/(1-p1))/(p0/(1-p0));
}ex16-2.stan
data {
int n;
int n0;
int n01;
int n1;
int n11;
}
parameters {
real<lower=0,upper=1> p0;
real<lower=0,upper=1> p1;
}
model {
n01~binomial(n0,p0);
n11~binomial(n1,p1);
}
generated quantities {
real RR;
RR=p1/p0;
real OR;
OR=(p1/(1-p1))/(p0/(1-p0));
real crv;
crv=(n0*(1-p0)*n1*p1-n0*p0*n1*(1-p1))/(n0*n1*(n0*(1-p0)+n1*(1-p1))*(n0*p0+n1*p1))^.5;
real k;
real po;
po=(n0*(1-p0)+n1*p1)/n;
real pe;
pe=n0/n*(n0*(1-p0)+n1*(1-p1))/n*(n0*p0+n1*p1)/n*n1/n;
k=(po-pe)/(1-pe);
}n0=30
n01=rbinom(1,n0,0.3)
n1=30
n11=rbinom(1,n1,0.6)
data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)
mdl=cmdstan_model('./ex16-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -48.2981
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -34.2445 0.00167309 0.000743161 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -34.24
## p0 0.17
## p1 0.47
## RR 2.80
## OR 4.38mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -38.62 -38.31 1.06 0.74 -40.79 -37.61 1.00 946 1146
## p0 0.19 0.18 0.07 0.07 0.09 0.31 1.00 1380 995
## p1 0.47 0.47 0.09 0.09 0.33 0.61 1.00 2037 1292
## RR 2.93 2.58 1.48 1.05 1.38 5.64 1.00 1423 1217
## OR 4.93 4.07 3.38 2.27 1.65 11.35 1.00 1487 1191
data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)
mdl=cmdstan_model('./ex16-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -52.7725
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 -34.2445 0.000159391 6.52948e-05 0.9447 0.9447 9
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -34.24
## p0 0.17
## p1 0.47
## RR 2.80
## OR 4.37
## crv 0.32
## k 0.63
## po 0.65
## pe 0.05mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -38.62 -38.31 1.06 0.74 -40.79 -37.61 1.00 946 1146
## p0 0.19 0.18 0.07 0.07 0.09 0.31 1.00 1380 995
## p1 0.47 0.47 0.09 0.09 0.33 0.61 1.00 2037 1292
## RR 2.93 2.58 1.48 1.05 1.38 5.64 1.00 1423 1217
## OR 4.93 4.07 3.38 2.27 1.65 11.35 1.00 1487 1191
## crv 0.30 0.31 0.12 0.11 0.12 0.49 1.00 1615 1415
## k 0.62 0.62 0.06 0.06 0.53 0.72 1.00 1671 1415
## po 0.64 0.64 0.06 0.06 0.55 0.73 1.00 1700 1404
## pe 0.05 0.05 0.00 0.00 0.05 0.06 1.00 1910 1637
point of subjective equality PSE
PSE: 50% threshold for sensing the difference between two stimuli is equal
JND: Just noticeable difference, difference between 50% threshold and 75%
r~B(n,p) #reaction for stimuli
p=1/(1+exp(-(a+b*x)))
x=x1-x0, x0,x1 is strength of stimuli
PSE=-a/b
JND=(log(0.75/0.25)-a)/b-PSEex6-3-0.stan
mulit logistic regression
data{
int N;
int m;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
}
model{
y~binomial_logit(m,b0+b1*x);
}n=30
m=10
x=runif(n,-2,2)
y=rbinom(n,m,1/(1+exp(-(-2+3*x))))
glm(y/m~x,family=binomial('logit'))##
## Call: glm(formula = y/m ~ x, family = binomial("logit"))
##
## Coefficients:
## (Intercept) x
## -1.79 2.84
##
## Degrees of Freedom: 29 Total (i.e. Null); 28 Residual
## Null Deviance: 22.5
## Residual Deviance: 2.72 AIC: 14data=list(N=n,m=m,x=x,y=y)
mdl=cmdstan_model('./ex6-3-0.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -694.414
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 8 -83.3571 0.00648888 0.00171875 1 1 10
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -83.36
## b0 -1.79
## b1 2.84mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -84.33 -84.05 0.95 0.68 -86.30 -83.42 1.01 784 850
## b0 -1.84 -1.82 0.29 0.29 -2.36 -1.39 1.00 630 845
## b1 2.92 2.89 0.36 0.36 2.34 3.57 1.01 612 870
b0=mcmc$draws('b0')
b1=mcmc$draws('b1')
pse=as.vector(-b0/b1)
quantile(pse,probs=c(0.0,0.025,0.05,0.25,0.5,0.75,0.95,0.975,1),na.rm=T)## 0% 2.5% 5% 25% 50% 75% 95% 97.5% 100%
## 0.426 0.508 0.522 0.587 0.630 0.676 0.741 0.762 0.861jnd=(log(0.75/0.25)-b0)/b1-pse
quantile(jnd,probs=c(0.0,0.025,0.05,0.25,0.5,0.75,0.95,0.975,1),na.rm=T)## 0% 2.5% 5% 25% 50% 75% 95% 97.5% 100%
## 0.262 0.297 0.308 0.349 0.380 0.413 0.469 0.485 0.562